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We study analytically the quantum thermalization of two coupled two-level systems (TLSs), which are con- 
nected with either two independent heat baths (IHBs) or a common heat bath (CHB). We understand the quantum 
thermalization in eigenstate and bare-state representations when the coupling between the two TLSs is stronger 
and weaker than the TLS-bath couplings, respectively. In the IHB case, we find that when the two IHBs have the 
same temperatures, the two coupled TLSs in eigenstate representation can be thermalized with the same tem- 
perature as those of the IHBs. However, in the case of two IHBs at different temperatures, just when the energy 
detuning between the two TLSs satisfies a special condition, the two coupled TLSs in eigenstate representation 
can be thermalized with an immediate temperature between those of the two IHBs. In bare-state representation, 
we find a counterintuitive phenomenon that, under some conditions, the temperature of the TLS connected with 
the high-temperature bath is lower than that of the other TLS, which is connected with the low-temperature bath. 
In the CHB case, the coupled TLSs in eigenstate representation can be thermalized with the same temperature as 
that of the CHB in nonresonant cases. In bare-state representation, the TLS with a larger energy separation can 
be thermalized to a thermal equilibrium with a lower temperature. In the resonant case, we find a phenomenon 
of anti-thermalization. We also study the steady-state entanglement between the two TLSs in both the IHB and 
CHB cases. 

PACS numbers: 03.65.Yz, 05.30.-d, 05.70.Ln 



I. INTRODUCTION 

Conventional quantum thermalization (ll |31 is understood 
as an irreversibly dynamic process under which a quantum 
system immersed in a heat bath approaches a thermal equi- 
librium state with the same temperature as that of the bath. 
The thermal equilibrium state |3i] of the thermalized sys- 
tem reads as pthiT) - Z^' txp{—Hs HkBT)}, which is merely 
determined by the Hamiltonian Hs of the thermalized sys- 
tem and the temperature T of its environment, where Zg - 
Trs {exp[-//s /(^sr)]) is the partition function of the thermal- 
ized system, with kg being the Bolztmann constant. During 
the course of a quantum thermalization, all of the initial in- 
formation of the thermalized system is totally erased by its 
environment. Recently, much attention has been paid to quan- 
tum thermalization (e.g., Refs. ll4l- [l2ll ). Specifically, a new 
kind of thermalization, called canonical thermalization (e.g., 
Refs. II13I - U7I1 ). has been proposed. 

The conventional quantum thermalization works for the sit- 
uations wherein one quantum system is connected with just 
one environment at thermal equilibrium. When we consider a 
composite quantum system, which is constructed with many 
subsystems and connected with many environments, the con- 
ventional quantum thermalization is no longer valid. There- 
fore, the density operator of the composite quantum system 
at thermal equilibrium can not be written as p,h{T). In fact, 
quantum thermalization of a composite quantum system is 
a very complex problem. On one hand, from the viewpoint 
of the environments, the composite quantum system could be 
connected with many independent environments or a common 



environment. At the same time, the temperatures of these en- 
vironments could be the same or different. On the other hand, 
the coupling strengths among these subsystems can affect the 
physical picture to describe the quantum thermalization of 
the composite quantum system. When the coupling strengths 
among these subsystems are stronger than the system-bath 
couplings, the composite quantum system can be regarded as 
a single system, while it is regarded as many individual sub- 
systems when the couplings among them are weaker than the 
system-bath couplings. 

The above mentioned situations come true in recent years 
since quantum systems can be manufactured to be more and 
more complicated and small, based on the great advances in 
physics, chemistry, and biology [18]. Therefore, the research 
on thermodynamics of small systems becomes very interest- 
ing. In particular, quantum thermalization of composite quan- 
tum systems becomes an important topic since many impor- 
tant results in this field, such as nonequilibrium work relations 
and fluctuation theorem (e.g., Refs. il9l - t23ll ). are based on the 
thermal equilibrium state of the composite quantum systems. 
As composite quantum systems are composed of many sub- 
systems, they could be connected either with many indepen- 
dent environments at different temperatures or a common en- 
vironment. Therefore, it is natural to ask the following ques- 
tions: (i) How do the couplings among the subsystems affect 
the quantum thermalization of a composite quantum system? 
(ii) What are the steady-state properties of a composite quan- 
tum system when it is connected with either many indepen- 
dent environments or a common environment? 

With these questions, in this paper, we study the quantum 
thermalization of two coupled two-level systems (TLSs) that 



are immersed in either two independent heat baths (IHBs) or a 
common heat bath (CHB). Simple as this model is, it is illus- 
trative. When the coupling between the two TLSs is stronger 
than the TLS-bath couplings, the two TLSs can be consid- 
ered as an effective composite system, i.e., a four-level sys- 
tem, and then we understand the quantum thermalization in 
eigenstate representation of the composite system. In addi- 
tion, when the coupling between the two TLSs is weaker than 
the TLS-bath couplings, we understand the quantum thermal- 
ization from the viewpoint of each individual TLS. However, 
due to the TLS-TLS coupling, the effective temperatures of 
the two TLSs should be different from those for the decou- 
pling case. 

As for the environments, there are two kinds of different 
situations: the IHB case and the CHB case. In the IHB case, 
we find that, when the two IHBs have the same tempera- 
tures, the two coupled TLSs in the eigenstate representation 
can be thermalized with the same temperature as those of the 
IHBs. However, in the case where the two IHBs have dif- 
ferent temperatures, just when the energy detuning between 
the two TLSs satisfies a special condition, the two coupled 
TLSs in eigenstate representation can be thermalized with an 
immediate temperature between those of the two IHBs. In 
the bare-state representation, we find a counterintuitive phe- 
nomenon that, under some conditions, the temperature of the 
TLS connected with the high-temperature heat bath is lower 
than that of the other TLS, which is connected with the low- 
temperature heat bath. 

In the CHB case, we also study the quantum thermalization 
in eigenstate and bare-state representations. In the eigenstate 
representation, the present case reduces to the conventional 
quantum thermalization, [i.e., one quantum system (an effec- 
tive four-level system formed by the two coupled TLSs) is 
thermalized by one environment (the common heat bath) in 
thermal equilibrium]. In addition, we also investigate the ef- 
fective temperatures of the two TLSs in bare-state representa- 
tion. It is found that the TLS with a larger energy separation 
can be thermalized with a lower temperature. In particular, 
in the resonant case, we find a quantum phenomenon of anti- 
thermalization when the two TLSs are connected with a com- 
mon heat bath. 

This paper is organized as follows: In Sec. |lll we present 
the physical models and their Hamiltonians. In Sec. Hill we 
study the quantum thermalization of the two coupled TLSs 
immersed in two IHBs. In Sec. IIVI we consider the case 
wherein the two TLSs are immersed in a CHB. We conclude 
this work in Sec.[Vl Finally, we give two appendices. lAlandlBl 
for detailed derivations of quantum master equations and tran- 
sition rates for the IHB and CHB cases, respectively. 



II. PHYSICAL MODELS AND HAMILTONIANS 

Let us start with introducing the physical models [as illus- 
trated in Figs.[na) and [Tib)]: two TLSs, denoted by TLSl 
and TLS2 with respective energy separations wj and W2, cou- 
ple with each other via a dipole-dipole interaction of strength 
^. At the same time, the two TLSs couple inevitably with the 
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FIG. 1: (Color online) Schematic diagram of the physical models 
under consideration. Two coupled two-level systems, denoted by 
TLSl and TLS2 with respective energy separations aii and ti;2, are 
connected with either (a) two IHBs at temperatures T[ and T2 or 
(b) a CHB at temperature T. Between the two TLSs, there exists 
a dipole-dipole interaction of strength ^. (c) The levels of the four 
bare states I;;,) (; = 1,2,3,4) of the two free TLSs. (d) The levels 
of the four eigenstates |/i,> (; = 1, 2, 3, 4) of Hamiltonian l[2j for the 
two coupled TLSs. In the presence of the baths, there exist bath- 
induced exciting and damping processes among the four eigenstates. 
The effective transition rates from states |/i,) to \Aj) are denoted by 
Tjj. Other cross-dephasing processes are not denoted explicitly. 



environments surrounding them. Specifically, in this paper, 
we consider two kinds of different cases: the IHB case and 
the CHB case. In the former case, the two TLSs are immersed 
in two IHBs, while, in the latter case, the two TLSs are im- 
mersed in a CHB. 

The Hamiltonian of the total system, including the two 
TLSs and their environments, is composed of three parts: 



H — H-T 



: + //«+ H, 



(1) 



where i/xLSs is the Hamiltonian of the two coupled TLSs, Hg 
is the Hamiltonian of the heat baths, and /// describes the in- 
teraction Hamiltonian between the two TLSs and their baths. 
The Hamiltonian //tlss (with Ti-Y) reads as 

Htlss = —(r\ + —(T^ + ^(o-|o-2 + 0-1 cr+). (2) 

The first two terms in Eq. (|2]) are free Hamiltonians of the two 
TLSs, which are described by the operators cr^ - {crj)' - 
\e)u{g\ and 07 = \e)ii{e\ - \g)u{g\ (l = 1,2), where \g)i and 
\e)i are, respectively, the ground and excited states of the /th 
TLS (i.e., TLSO- The last term in Eq. (|2| describes a dipole- 
dipole interaction of strength ^ between the two TLSs. We 
note that the Hamiltonian given in Eq. ^ has been widely 



studied in various physical problems, such as quantum logic 
gates t24.1 . coherent excitation energy transfer p3\, decoher- 
ence dynamics ll26l - l28ll . and nonequilibrium thermal entan- 
glement HI, [Ml ■ 

The Hilbert space of the two TLSs may be spanned by the 
following four bare states I771) = \ee), \t]2) = \eg), It/s) = \ge), 
and I774) = 1^^) [as shown in Fig.fTfc)], which are eigenstates 
of the free Hamiltonian (wicr-j + (y20^)/2 of the two TLSs, 
with the corresponding eigenenergies £■,„ = -E,j^ = (d,„ and 
E,p = -Enj = Aw/2. Here we have introduced the mean 



the two TLSs are immersed in a CHB with the Hamiltonian 

^(CHB) ^ ^ oj^ja^.aj, (6) 

J 

where a and Oj are, respectively, the creation and annihilation 
operators of the jth harmonic oscillator with frequency ojaj. 
The interaction Hamiltonian between the two TLSs and the 
CHB reads as 



(wi + W2)/2 and the energy detuning //f ™' = ^ 8ij(a\a-^ + a-\aj) + ^ gij{a\cr^ + (rlaj). (7) 



energy separation u, 

Aw — U)l - U)2. 

Due to the dipole-dipole interaction, a stationary single- 
excitation state should be delocalized and composed of a com- 
bination of the single-excitation states in the two TLSs. Ac- 
cording to Hamiltonian (|2]i, we can solve the eigenequation 
//tlssI'Im) = E,i,}A„), (n - 1,2,3,4) to obtain the following 
four eigenstates [as shown in Fig. [lid)]: 

|/li ) = 1/71 >, 1^2) = cos(0/2)|772> + sin(0/2)|773), 
M3) = - sin(0/2)|)72> + cos(0/2)|773), 1^4) = \m)^ (3) 

with the corresponding eigenenergies E^, - -E,i^ = w„, and 
£^1, - -£^3 - s/AaP-IA + ^. Here we have introduced a 
mixing angle 6 {0 < 6 < n) hy tan0 = 2^/ Aw. For a 
positive ^, when wi > W2, namely Aw > 0, we choose 
- arctan(2^/Aw); When wi < W2, that is Aw < 0, we 
choose 9 - n + arctan(2^/Aw). The dipole-dipole interaction 
mixes the two bare states I772) and I773) with one excitation, and 
does not change the bare states I771) with two excitations and 
I774) with zero excitation. 

Aside from the dipole-dipole interaction between the two 
TLSs, there exist couplings between the TLSs and their en- 
vironments. In general, when the couplings of a system with 
its environment are weak, it is universal to model the environ- 
ment as a harmonic-oscillator heat bath and choose a linear 
coupling between the system and its environment [31]. In this 
paper, we consider this situation and model the environments 
as harmonic-oscillator heat baths. As mentioned above, we 
will consider two kinds of different cases: one is the IHB case 
and the other is the CHB case. 

In the IHB case, as shown in Fig. [Ha), the two TLSs are 
immersed in two IHBs described by the Hamiltonian 



H' 



(IHB) 



^i/i^'+i/' 



ib) 



H' 



(a) 



-!■ 



Jaf 



a -Oj 

J ^ 



H' 



ih) 



^Yj^bkblbk- (4) 



iW 



f(b) 



Here H^ and H^ are, respectively, the Hamiltonians of the 
baths for TLS 1 and TLS2. The creation and annihilation oper- 
ators a . (bj^) and Oj (bk) describe the jth (kth) harmonic oscil- 
lator with frequency coaj icohk)- The interaction Hamiltonian 
of the two TLSs with the two IHBs reads 



^/™^ = Yj SiA^y\ + ^l«i) + Yj Sikibla-^ + cr^bk). (5) 



On the other hand, in the CHB case, as shown in Fig.lTJb), 



For simplicity, we have assumed that the TLS-bath coupling 
strengths gij, g2j, and g2k are real numbers. 



III. QUANTUM THERMALIZATION OF TWO COUPLED 
TLSS IMMERSED IN TWO IHBS 

In this section, we study the quantum thermalization of the 
two coupled TLSs that are immersed in two IHBs. We depict 
the evolution of the two TLSs in terms of a quantum master 
equation. By solving the equations of motion of the density 
matrix elements to obtain steady-state solution, we study the 
steady-state properties of the two coupled TLSs. 



A. Equations of motion and steady-state solutions 

We consider the situation wherein the environments of the 
two TLSs are memory-less and the couplings between the 
TLSs and the environments are weak. Then, we may adopt 
the usual Born-Markov approximation in derivation of quan- 
tum master equation. At the same time, we derive the mas- 
ter equation in eigenstate representation of the two coupled 
TLSs so that we can safely make the secular approximation 
(equivalent to the rotating-wave approximation) to obtain a 
time-independent quantum master equation by neglecting the 
rapidly oscillating terms |1]. Therefore, our discussions are 
under the Born-Markov framework. 

In the case of two IHBs, the evolution of the two coupled 
TLSs is governed by the following Born-Markov master equa- 
tion in the Schrodinger picture, 

Ps = i[ps , H-Yi^ss] 

+ X [rjv(2T,7Psr^-; - Tjjps - PsTjj) 

+Tij{2TjiPsTij - TuPs - PsTii)\ 

-I-2Ai(T42PsTi3 + T31PST24) 

+2A2(T2iPsT34 + T43PST12) 

+2A3(T24PsT31 +T13PST42) 

-I-2A4(ti2PsT43 -I-T34PST21), (8) 

which wiU be derived in detail in AppendixlAl In Eq. ([8]), ps is 
the reduced density operator of the two TLSs. The operators 
Tij are defined by r^ - \Ai){Aj\ in the eigenstate representation 



of the two coupled TLSs. The summation parameters (/, j) in 
the second line of Eq. ^ can take (/, ;) = (4, 2), (3, 1), (2, 1), 
and (4, 3). In the present model, the effective rates in Eq. (O 
are defined as Fn = r24 = Ei, r3i = F42 = E2, E12 = F34 - E3, 
and E21 - E43 = E4, with 



El 


= cos2(6»/2)Ai(ei) + sm\6/2)Bi(ei), 


r2 


= cos2(6»/2)A2(ei) + sm^ie/2)B2isi), 


r3 


= sin2(6»/2)Ai(e2) + cos^{0/2)Bi(e2), 


E4 


= sin2(6»/2)A2(e2) + cos2(6i/2)B2(e2), 


Ai 


= cos^(6'/2)Ai(ei) - sin^(0/2)Bi(ei), 


A2 


= - sin2(0/2)Ai(e2) + cos2(6'/2)Bi(e2) 


A3 


= cos2(6'/2)A2(£i) - sin2(6i/2)B2(ei), 


A4 


= - sin2(0/2)A2(e2) + cos2(6'/2)B2(e2) 



(9) 



where we define Ai(e,) - yaCe,) [«„(£,) + 1], A2(e,) - 
7a(ei)naiEi), Bi(si) = jbisdlnbisi) + 1], and B2(e,) = 
yb{Ei)hh{ei) (i = 1,2). The energy separations ei and S2 are 
introduced as si = £,1, --E^, - -E^, -EAi - (^m + yjAcdP-fA- + ^ 
and 62 = -£■,!, - -E^, - £^, - -E^^ - (jj„, - -\JAcj^/4- + ^^. The 
rates y„(e/) = ^0a(si)g](si) and jhisi) = 7Tgb(Ei)gl{Ei), where 
Qaisi) and Qbisi) are, respectively, the densities of state at en- 
ergy e, of the heat baths for TES 1 and TES2. In the following 
we assume ya(si) = yb(si) = n and yaisi) = yhisi) = Ta- 
in addition, we introduce the average thermal excitation num- 
bers «fl(e,) - l/[exp{si/Ti) - 1] (with the Boltzmann constant 
ks = 1) and ni(e,) - l/[exp{si/T2) - 1] (/ = 1, 2) for the heat 
baths of TES 1 and TES2, at the respective temperatures Ti 
andralUIl. 

Based on quantum master equation ([8]l, it is straightforward 
to obtain optical Bloch equations for the density matrix ele- 
ments of the two TESs in the eigenstate representation. De- 
noting X(0 = «Tii(0), <T22(f)), <T33(f)), <T44(0))^ CT" Stands 
for matrix transpose), then the optical Bloch equations for the 
diagonal density matrix elements in the eigenstate representa- 
tion can be expressed as 



X(0 = M^™^X(0, 



where the coefficient matrix M'™^^ is defined by 



M^HB) = -2 



(10) 



( El + E3 -E4 


-E2 


-E3 El + E4 


-E2 


-E, 


E2 + E3 -E4 


I -El 


-E3 E2 + E4 



(11) 



From Eq. ( fTOl i. we can see that the evolution of the diago- 
nal density matrix elements decouples with off'-diagonal ele- 
ments. 

The transient solutions of optical Bloch equation (fTOl l can 
be obtained with the Eaplace transform method. To study 
quantum thermalization, however, it is sufficient to obtain the 



steady-state solutions, 

{T22} ss - 

<T33)ss - 

{T44) ss - 



E2F 



2^4 



(Fi+F2)(F3+F4)' 

r2F3 

(Ei+F2)(E3+F4)' 

riF4 

(Ei+F2)(E3+F4)' 

E1E3 

(Fi+F2)(F3+F4)' 



(12) 



where the subscript "ss" means steady-state solutions. 

According to Eq. (jS), we can also obtain the equations of 
motion for these off-diagonal density matrix elements of the 
two TESs as follows: 

<-r2i(0) = -(2Fi+F3+F4-/e2)<r2i(0>+2A3<r43(f)), 

<i-3i(0> = -(Fi + E2 + 2E3 - /ei)<r3i(0) + 2A4<T42(0), 

<T-4i(0) = -(El + E2 + E3 + F4 - isi - ie2){T4i(t)}, 

(fn(t)) = -(Fi + F2 + E3 + F4 - isi + /e2)<T32(f)), 

<t42(0) = -(Fl + F2 + 2E4 - iei){T42(t)) + 2A2<T3l(f)), 
<f43(f)) = -(2F2+F3+F4-!e2)<T43(0) + 2Ai<T2l(0>. 

(13) 

Other off-diagonal elements can be obtained via (Ty(f))* = 
{Tji{t)}. It can be found that the steady-state solutions of the 
equations of motion for these off-diagonal elements are zero. 



{Tij)s.9 = 0, / ^ ;'. 



(14) 



Based on the steady-state solutions for these density matrix 
elements, we can analyze the steady-state properties of the 
two coupled TESs. 



B. Quantum thermalization in eigenstate representation 

For the present system with two IHBs, when the coupling 
between the two TESs is stronger than the TES-bath cou- 
plings, the two coupled TESs can be regarded as an effec- 
tive four-level system connected with two IHBs. Therefore, 
in the eigenstate representation, the dynamic evolution pro- 
cess of the two-coupled TESs approaching their steady state of 
thermal equilibrium can be understood as a non-equilibrium 
quantum thermalization: thermalization of a quanturn system 
connected with many IHBs at different temperatures i32l - l39ll . 

As the steady state of the two coupled TESs is completely 
mixed in eigenstate representation, we can introduce effec- 
tive temperatures to characterize the relation between any two 
eigenstates based on their steady-state populations. From 
Eq. (fT2b we can find the following relations 



(15) 
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FIG. 2: (Color online) Plot of the scaled effective temperatures 
Tefi(£i)/~y (solid blue line) and T^ff(s2)/y (dashed red line) vs the 
mixing angle 9. Other parameters are set as 71 = 72 = 7, f/y = 10, 
w,„/7 = 20, TJy = 5, and Ti/y = 10. 



Generally, it is impossible to define an effective temperature 
for the effective four-level system at steady state. Accord- 
ing to Eq. (flSl l. we can characterize the state of the effective 
four-level system via introducing the following two effective 
temperatures ll40ll by 



<Tll)., 
<T22)i 



(T33)s.; 
{t 44) ss 

= exp 



exp 



7'eff(e2) 



£1 -£2 



<T"33> 

In terms of Eq. ( fTSl l. we get 



Teffisi - ei) 



Tes{e2) 



in(r,/r4)- 



(16) 



(17) 



When the two IHBs have the same temperatures, namely 
T\ - T2 - T, we can show that the above two effective tem- 
peratures are equal to those of the two IHBs, that is 



7'eff(ei) = 7'eff(e2) = T. 



(18) 



For the nonequilibrium case of Ti + To, the two effective 
temperatures are different for general case. We find the fol- 
lowing relations min(ri,r2) < T^s{e\) < max(ri,r2) and 
min(ri,r2) < retf(e2) < max(ri,r2), which mean that the 
two effective temperatures reff(£i) and reff(s2) will be within 
the region from min(ri, T2) to max(ri, T2) lUlll . Under some 
special conditions, the two effective temperatures could be 
equal. In this case, we consider that the effective four-level 
system is thermalized by the non-equilibrium environments: 
two IHBs at different temperatures. 

In Fig. |2] we plot the two effective temperatures reff(ei) 
and reff(e2) as a function of the mixing angle 6. It shows 
that reff(ei) and reff(e2) can be equal for a special mixing an- 
gle 00, which is determined by 



in(r,/r2) ~ in(rj7r4)' That is to 
say the effective four-level system formed by the two coupled 
TLSs can be thermalized with an effective temperature when 



Aw = ^/ tan 00. Note that this effective temperature actually 
is a nonequilibrium temperature hHi . We can also see from 
Fig.[2]that the two effective temperatures reff(ei) and Te«ie2) 
are within the region from Ti to T2. We point out that the 
mixing angle should be chosen to make sure that £1, > E^^. 



C. Quantum thermalization In bare-state representation 

When the coupling between the two TLSs is weaker than 
the TLS-bath couplings, we understand the quantum thermal- 
ization in bare-state representation. We can express the bare 
states with the eigenstates as 

1^1) = 1^1), I772) = cos(0/2)|^2) - sin(0/2)|^3), 

I773) = sin(0/2)|i2) + cos(0/2)|i3), \ri4)^\M)- (19) 

Then we can obtain the relations 

<f^^l_2(0) = <Tii(0)-(T44(0) 

-(-1)'COS0[<T22(O)-<T33(O)] 

+ (-l)'sin0[<T23(f)) + <T32(O)]. (20) 

According to Eqs. ( fT2l l. ( fT4l l. and (l20t . the steady-state solu- 
tion can be obtained as 



<oi_^, 2)"" = 



(r2r4-rir3)-(-i)'cos0(r2r3-rir4) 



(Fi + r2)(F3 + F4) 



(21) 



In addition, the off-diagonal elements of the density matrices 
of the two TLSs can be expressed as 






Sin(0/2)[<T14(O) - <T34(0)] 
+ COS(0/2)[<Ti3(O) + <T24(O)], 
Sin(0/2)[<T24(O) - <T13(0)] 
+ COS(0/2)[<Ti2(O) + <T34(O)]. 



(22) 



Because of {Tij(t))ss -Q{i + j), we have {cr^)ss - (C2 ^.s.s - 0, 
which implies that the steady states of the two TLSs in bare- 
state representation are completely mixed. Based on these, it 
is possible to introduce two effective temperatures of the two 
TLSs as follows: 



7'eff(w/) = 



CO, 






/= 1,2. 



(23) 



In 



In Fig. [3] we plot the effective temperatures TegicJi) and 
Tef{{a>2) as a function of the mixing angle 9 for various tem- 
perature distributions. From Fig. [3] we can see the following 
three interesting results. Firstly, when 6 x n/l, Teff(cL)\) and 
Tes{a>2) become approximately equal. Specially, at the reso- 
nant point 9 = 7r/2, the two TLSs have the same temperatures 
no matter whether the two bath temperatures are the same or 
not (see the cross point at 9 = n/2 in figures). This result 
can be explained as follows: When the two TLSs are nearly 
in resonance with each other, the dipole-dipole interaction can 




FIG. 3: (Color online) Plot of the scaled effective temperatures 
Tcff{(jjC)ly (solid blue line) and T^ff{(jj2)ly (dashed red line) vs the 
mixing angle 6 for various temperature distributions (shown in fig- 
ures). Other parameters are set as 71 = 72 = 7. ^/7 = 0.1, and 

oj,„ly = 20. 



induce population exchange between the two TLSs such that 
their temperatures are approximately equal. When = n/l 
{oji = 0)2), Eq. (|2TJ reduces to 



{o-\)ss - {0-7} ss 



r2r4-rir3 

(ri+r2)(r3 + r4)' 



then we have Tef{(aji) = Teff((L)2). 

Secondly, when 6 is near to and n, reif(wi) and reif(w2) 
approach approximately Ti and T2, respectively. This result 
can be understood as follows: Near to 6 x Q and n, the en- 
ergy detuning \A(l)\ - 2^/\ tan d\ between the two TLSs is very 
large, and the population exchange between them can be ne- 
glected for a weak coupling strength ^. At the same time, the 
energy separation shifts of the two TLSs are neglectable [this 
can be seen from the following effective Hamiltonian (IZSlll. 
Hence, the temperatures of the two TLSs should be equal to 
those of their IHBs, respectively. 

Thirdly, in Fig. [3] there are some regions of 6 where the ef- 
fective temperature Tf^fficoi) of TLS 1 could be smaller than the 
temperature Tsf{(a)2) of TLS2 although the bath temperatures 
Ti > T2- This is a counterintuitive result. Intuitively, in these 
cases Tfifi((L)2) should be smaller than 7'eff(wi), because TLS2 
is connected with a low-temperature bath while TLS 1 is con- 
nected with a high-temperature bath [In Fig. [3ja), when the 
two bath temperatures are the same, the effective temperatures 
of the two TLSs should be equal]. Actually, the counterintu- 
itive result is a net effect of the energy separations of the two 
TLSs, the bath temperatures, and the coupling between the 
two TLSs (coupling induces population exchange and energy 
shifts). From Figs. |3la) to|3jd), the curves change gradually 
with the increase of the temperature difference T\-T2 between 
the two baths. And the counterintuitive region decreases with 
the increase of the temperature difference Ti - 7^2. 

In the following, we present a microscopic explanation 
for Fig. [3] as a physical insight of the counterintuitive phe- 
nomenon. When the values of 9 are far from the resonant 



point (not near to and n), the two TLSs will have large 
detuning from each other. Under the large detuning con- 
dition ^ <s A(jj, the real population exchange between the 
two TLSs is compressed, the dipole-dipole interaction induces 
energy shift to the two TLSs. In this case, we can derive 
an effective Hamiltonian to describe the two TLSs with the 
Frohlich-Nakajima transformation approach 11421 14311 . Start- 
ing from the Hamiltonian //jlss - H^ + H'j with H'^ - 
{(i)\(T\ + W20^)/2 and H'j - ^(cr|cr^ + cr^crl^), we introduce an 
operator S - -^{o-^cr:^ ~ cr'^cr^) j Au, which meets the condi- 
tion H'j + [H'q,S] - 0. Then the effective Hamiltonian reads 
as 



//eff = //o + -[H'i,S] = 2^\cr'i + 2^^20-1' 



(25) 



where the shifted energy separations are defined hy (b\ - (jJ\ + 
^ lAui and Uj2-<j^2- ^ I Aoj. 

We can see from Eq. ( |25] | that, under the large detun- 
ing condition, there is no effective coupling between the two 
TLSs. The dipoIe-dipoIe interaction between the two TLSs 
shifts their energy separations slightly. Hence, when the two 
TLSs (with the shifted energy separation) are thermalized to 
thermal equilibrium with their baths, we have the relation 
exp(-w,/r,) = pflp-^ (^ = 1,2) for the TLS/ (with shifted 
energy separation w;, excited and ground state populations p\ 
and p^) in thermal equilibrium at temperature T/, and then the 



(24) effective temperatures defined in Eq. ( |23T l should be 



7'eff(w;) = 



W, 






Upflpf) "^i 



(26) 



From Eq. (l25T l. we can see that, for a positive ^, when < < 
7r/2, we have Auj > 0, then wi > wj and d>2 < (x>2- Hence 
the effective temperatures T^fficoi) < T^ and T^ff{(jj2) > T2- 
On the other hand, when 7r/2 < 6 < n, we have Aw < 0, 
then (bi < a>i and 0*2 > C02- Hence, the effective temperatures 
TefficJi) > Ti and Teff{cJ2) < 7^2. 

According to the above analysis, we can see that, when 
T\ - T2 [Fig.[3ja)], there will exist a counterintuitive region. 
At the same time, when 6 is near to and n, the shifted energy 
separations^ /Aw approaches zero, then Tsff(cji) « T/. Hence, 
with the increase of the bath temperature difference Ti - T2, 
the difference between the two effective temperatures also in- 
crease, which leads to the counterintuitive region decreases. 
These results can be seen from Fig. [3] 

In fact, the above intuitive result is based on the phe- 
nomenological master equation 



with 



PS = i[ps, -H^TLSs] + -Clips] + -C2[Ps], 



71, 



(27) 



-Ci=i,2[ps] = y(«; + l)(2cr, pcr+ - o-fcr, p - pcrlo-, ) 

+—hi(2o-lpcrJ - o-JctIp - pcrjcrf). (28) 

The superoperator X.i[ps] describes the dissipation of a TLS/ 
(/ - 1,2) immersed in a heat bath at temperature T; (hi = 
l/[exp(cji/Ti) - 1]). Therefore, Eq. ( l27T i is not valid in the 



case of two coupled TLSs, especially when the coupling be- 
tween the two TLSs is stronger than the TLS-bath couplings. 
In addition, in Eq. (l27T i. the effects on the TLSs from the two 
baths are different, one is direct and the other is indirect. For 
example, the bath of TLS 1 affects TLS 1 directly, while the 
bath of TLS2 affects TLSl indirectly through TLS2. On the 
contrary, our results given in Eq. ( |23] ) are based on quan- 
tum master equation dS), which is rigorously derived in the 
eigen-representation of the two coupled TLSs. Hence, the 
dissipation is depicted in the eigen-representation of the two 
coupled TLSs. In other words, the two TLSs play equiva- 
lent roles and the two baths directly affect the TLSs. The 
resonant case is a clear example for the equivalent role of 
the two TLSs. In the resonant case 6 - njl, we obtain 
Tefi{<jJ\) - Teg((jj2) in terms of quantum master equation (O, 
while we get Tsf{(a>i) > Tsff(ciJ2) from Eq. (l27T i when Ti > T2- 



D. Steady-state entanglement between the two TLSs 

In the IHB case, there exists a dipole-dipole interaction be- 
tween the two TLSs. Therefore, a natural question is: what 
the quantum entanglement is between the two TLSs after they 
are thermalized. As we know, during the thermalization pro- 
cesses (not anti thermalization), all of the initial information of 
the two coupled TLSs is totally erased, and the steady state of 
the two TLSs is determined by the decay rates and bath tem- 
peratures. Hence, we need to know the steady-state entangle- 
ment in the two TLSs. We note that entanglement dynamics in 
similar systems has been studied 1441 . In the following we ap- 
ply the concurrence to quantify the steady-state entanglement 
in the two TLSs. 

For a 2 X 2 quantum system (two TLSs) with a density ma- 
trix p expressed in the bare-state representation, its concur- 
rence [45] is defined as 



C(p) = maxjO, V^ " V^ ~ V^ ~ V^)' 



(29) 



where Sj (/ = 1,2,3,4) are the eigenvalues (si being the largest 
one) of the matrix pp. The operator p is defined as 

(30) 



p = (cr\ ® crl)p*(cr\ (g) cr^). 



where p* is the complex conjugate of p and cr^i is the Pauli 
matrix of TLS/. For the 2x2 system, C = and C - I mean, 
respectively, the density matrix p is an unentangled state and 
a maximally entangled state. In particular, for the so-called 
"X"-class state with the density matrix (expressed in the bare- 
state representation) 



(31) 
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the concurrence is W 



C(p) = 2max{0, |p23l - VPiiP44, \pu\ - VP22P33} ■ (32) 
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FIG. 4: (Color online) Plot of the steady-state concurrence C(p„) 
vs the scaled bath temperature T/y for various values of ^/y. Other 
parameters are set as 71 = 72 = 7, = ?r/2, Ti = T2 = T, and 

(^m/y = 20. 



Now, for the present system, its density matrix elements 
in bare-state representation can be expressed as {rij\p\T]i) - 
Tr[|77,)(77j|p] = Tr[jt/,jp] = {^ij) with the transition operator 
/^O - VldiiljV The density matrix elements in the eigenstate 
representation are expressed by {Aj\p\Ai) - Tr[|/l,)(/lj|p] = 
Tr[r,jp] = {Tij(t)} with t,j = |/l,)(/lj|. Since the concurrence 
is defined in the bare-state representation, and the evolution 
of the system is expressed in the eigenstate representation. 
Therefore we need to obtain the transformation between the 
two representations as follows: 



{finit)} = <Tii(f)), <A'44(f)> = <T44(0), 

{^I22(t)) = cos2(0/2)<r22(O> + sin2(0/2)<T33(f)>, 

--sin0[<T23(O) + <^32(O)], 

</l33(0) = Sm\0/2){T22(t)) + COs2(0/2)<T33(f)), 

+ isin0[(T23(/)) + (T32(O)], (33) 



and 



</Z23(0) = -sin2(0/2)<T32(O) + COs2(0/2)<T23(f)), 
+ -sin0[<T22(O)-<T33(O)], 

<//i2(0) = cos(0/2)<Ti2(f)) - sin(0/2)<ri3(f)>, 
<//i3(0> = sin(0/2)(ri2(O) + cos(0/2){t nit)), 

<jUi4(0) = <T14(0), 

(//24(f)> = COS(0/2)<T24(O> - Sin(0/2)<T34(O), 

<jU34(0) = Sin(0/2)<T24(O) + COS(0/2)<T34(f)>. (34) 

According to the steady-state solutions given in Eqs. ( fT2b 
and ( fT4b . the steady-state density matrix of the two TLSs in 
the bare-state representation can be obtained with the follow- 



ing nonzero elements 

(/Ull).s.! = {ti\)ss, {H44)ss - {T44>j.s, 

{m)ss = COs2(0/2)<T22).. + Sin2(6l/2)<T33). 
{^I33)ss = Sin2(0/2)<T22).. + COS2(0/2)<T33). 



1 



<y"32).« = 2 ^^"^^^^22)s 



The concurrence of this steady state is 

C(pss) = 2 max {o, K/"32)sj| - ^|{^inhA^M4} 



<T33).,.). (35) 



(36) 



In the following, we study the steady-state concurrence of 
the two coupled TLSs in the case of = n/2 and T\ - T2 - 
T. In Fig. m we plot the steady-state concurrence C(pss) as 
a function of the bath temperature T for various values of the 
dipole-dipole interaction strength ^. Figure |4] shows that a 
larger steady-state concurrence can be created for a larger ^. 
We also see the sudden death of the concurrence when the 
bath temperature increases up to a critical value. Notice that 
the phenomenon of threshold temperature has also been found 
in thermal entanglement of spin model ll29ll30ll47ll48ll . 



IV. QUANTUM THERMALIZATION OF TWO COUPLED 
TLSS IMMERSED IN A CHB 

In the above section, we have studied the quantum thermal- 
ization of two coupled TLSs immersed in two IHBs. However, 
in some cases, the two coupled TLSs can be considered to be 
placed in a CHB. In this section, we shall study the quantum 
thermalization of two coupled TLSs immersed in a CHB. 



A. Equations of motion and steady-state solutions 

The quantum master equation describing the evolution of 
the two TLSs immersed in a CHB at temperature T has the 
same form as Eq. ^, but the effective rates are not the same 
as those in the IHB case. In the CHB case, the rates read as 



ri2 = 

r2i = 

ri3 = 

r3i = 

r24 = 

Ui = 

r34 = 

r43 = 

Ai = 

A2 = 

A3 = 

A4 = 



sin(0/2) Vri(e2) + cos(0/2) ^fi^^] [n(e2) + 1], 
sin(0/2) Vri(e2) + cos(0/2) ^fM^f Hsi), 
cos(0/2) Vri(ei) - sin(0/2) yfM^f [«(ei) + 1], 
cos(0/2) Vri(ei) - sin(0/2) yf^^f n(ei), 
cos(6»/2) Vri(ei) + sin(6l/2) Vr2(ei)]^ [n(ei) + 1], 
cos(0/2) Vri(ei) + sin(0/2) ^fM^f h(si), 
sin(0/2) Vnfe) - cos(0/2) yf^^f [n(s2) + 1], 
sin(6i/2) Vri(e2) - cos(0/2) ^Jy2(s2)f n(e2), 

cos2(0/2)ri(ei) - sin2(0/2)r2(ei)][«(ei) + 1], 

• sm\0/2)yi(s2) + cos\e/2)y2(e2)][h(s2) +1], 
cos\0/2)yi{si) - sm\0/2)y2(ei)]n(ei), 
- sm\e/2)yi(e2) + cos\6/2)y2{s2)]n(e2), (37) 



where we introduce the rates yiisj) - 7Tg{si)g^(si), y2{si) = 
nQ{Ei)g\{E{), and the average thermal excitation number 
n(e,) — l/[exp(e,7r)- 1] (/ = 1,2). The detailed derivation of 
these rates will be given in Appendix IBJ In comparison with 
the rates in the IHB case, the rates in the CHB case have cor- 
relation terms which are induced by the CHB. In the following 
discussions, we assume yi(e,) = yoi^i) = y(si)- 

Correspondingly, the optical Bloch equation for the CHB 
case has the same form as Eq. ( flOb , but the coefficient matrix 
is replaced by M**--^^^ with the following expression 



M^CHB) 



( ri2 + To -r2i -r3i o 

-ri2 r2i + ^24 -r42 

-fn r3i + r34 -r43 

-r24 -r34 r42 + r43 j 

(38) 
Similar to the IHB case, the equations of motion for diagonal 
density matrix elements decouple with the off-diagonal ele- 
ments. The steady-state solutions of the present optical Bloch 
equation read as 



(t11 ).s.s - 

{^12)55 - 

{Ti-i)ss - 

(T44 )ss - 



(r2i + r24)r3ir43 + (r3i + r34)r2ir42 



(ri2 -1- ri3)r34r42 + (r42 + r43)ri2r3i 

A 

(ri2 + ri3)r43r24 + (r42 + r43)r2iri3 



(r2i + r24)ri3r34 + (Tj,\ + r34)ri2r24 



(39) 



with A = (ri2 + ri3)(r34r42 + r43r24) + (r2i + r24)(r3ir43 -i- 
ri3r34) -I- (r3i + r34)(r2ir42 + ri2r24) + (r42 + r43)(r2iri3 -i- 
ri2r3i). 

We can also obtain the equations of motion for these off- 
diagonal density matrix elements as follows: 

<t2i(0) = -(r2i+ri2 + ri3 + r24-/e2)<T2i(0) 

-h2A3<T43(0), 

<t3i(0) = -(ri2 + r3i+ri3 + r34-i£i)<T3i(f)> 

-h2A4(T42(0), 

<t4i(0) = -(ri2 + Fn + r42 + r43 - /ei - /e2)<T4i(0), 

(t32(0) = -(r21 + F31 + F24 + F34 - isx + ie2){T32(t)), 
<t42(0) = -(r21 + F42 + F24 + F43 - /ei)<T42(0) 

+2A2<T3i(0>, 
<t43(0) = -(F31 + F42 + F43 + F34 - /e2)(T43(0) 

+2Ai<T2i(r)). (40) 

The equations of motion for other elements can be obtained 
by {Tijif)) - {T*..(t)). Obviously, the steady-state solutions of 
these off-diagonal density matrix elements are zero. 



B. Quantum thermalization in eigenstate representation 

Differently from the IHB case, for the present four-level 
system, we introduce three effective temperatures to charac- 
terize its state. The three effective temperatures T\2, Tu, and 



T-i4 are defined according to the populations of the four levels 
as follows: 



10.002 



Tij = 



in(lfe) 



(41) 



We can show that the above introduced temperatures are the 
same as that of the CHB, 



Ti2 = Tn = 7^34 = T, 



(42) 



which means the two coupled TLSs can approach a thermal 
equilibrium with the CHB. In other words, the two coupled 
TLSs in eigenstate representation can be thermalized by the 
CHB. 

According to Eqs. (fTST l and ( l42l i. we know that when the 
temperatures of the heat baths are T, irrespective of two IHBs 
or a CHB, the effective four-level system formed by the two 
coupled TLSs can be thermalized into a thermal equilibrium 
state with the same temperature T . In other words, based on 
the thermal equilibrium state at temperature T, we can not 
know whether the two coupled TLSs are connected with two 
IHBs or a CHB. 



C. Quantum thermalization in bare-state representation 

We also investigate the quantum thermalization of the two 
TLSs in the bare-state representation. In terms of Eqs. (|20] | 
and ( [39] l. we can obtain the steady-state average values of the 
two Pauli operators cr^ and cr^ as follows: 



<0l=l,2>" = 



(Fai + r24)(r3ir43 - ri3r34) 



(r3i + r34)(r2ir42 - ri2r24) 
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FIG. 5: (Color online) Plot of the scaled effective temperatures 
Tcff(<^[)/7 (solid blue line) and T^ff(a)2)/'y (dashed red line) vs the 
mixing angle 6. Other parameters are set as -/(si) = y(s2) = y, 
^/y = 0.1, (i>„,/y = 20, and T/y = 10. The resonant point 6 = n/2 is 
useless. 



D. Quantum anti-thermalization in the resonant case 



When = n/2, the two TLSs are in resonance, and then 
the decay rates Fn, F31, F34, and F43 are zero under the as- 
sumption jiisi) = y2(ei)- Hence, the eigenstate \Aj,} decou- 
ples with other eigenstates, resulting in an anti-thermalization 
phenomenon. The state I/I3) is called as a "dark state" |49[]. 
In this case, we need to rewrite new optical Bloch equations 
for these density matrix elements. We obtain the equations of 
motion for diagonal density matrix elements 



(F12 + Fi3)(r34F42 - r43r24) 



-H(-1)'"'C0S( 



(F42 + F43)(Fi2F3i - F21F13) 



(43) 



Moreover, we have {o-'^)ss = (cj^s.? = 0. Similar to Eq. 
in the above section, we also introduce two effective temper- 
atures to characterize the state of the two TLSs. In Fig.|5] we 
plot the two effective temperatures as a function of the mix- 
ing angle 6. We emphasize that the resonant point = n/l in 
Fig.|5]should be taken out. 

We can draw a conclusion from Fig.|5]that the TLS with a 
larger energy separation can be thermalized to a thermal equi- 
librium state with a lower temperature. This can be seen as 
follows: When wi > 1^2, we have < < jt/2, from Fig.|5]it is 
clear that the temperature of TLSl is lower than that of TLS2. 
On the other hand, when a)[ < 0)2, we have 7r/2 < < n, 
Fig.|5]indicates that the temperature of TLS2 is lower than that 
of TLSl in this region. The physical explanation for Fig.|5]is 
the same as that for Fig.[3ja). 



<f„(0) = -2Fi2<Tn(0)+2F2i<T22(0), 

<f22(0) = 2Fi2<T„(0)-2(F2i+F24)<r22(0) 

+2F42<T44(0) 
<t33(0) = 0, 
<t44(0) = 2F24<T22(0) - 2F42<T44(0), (44) 



and for off-diagonal density matrix elements 



(T2l(t)} = -(F21 + F12 + F24 - /e2)<T2l(0), 

<t3l(0) = -(ri2-/ei)(T3l(0), 

<t4i(0) = -(F12 + F42 - iei - ie2){T4i(t)}, 

<T"32(0) = -(r21 + F24 - iSl + iS2){Ti2{t)), 

<t42(?)) = -(r21 + F42 + F24 - /ei)<T42(0), 

<f43(0) = -(r42 - /e2)<T43(0)- 



(45) 
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The steady-state solutions for these density matrix elements 



are 



{T22}ss 
{T44)ss 



[l-<T33(0))]r2ir42 

ri2r42 + ri2r42 + r2ir42 

[l-<T33(0))]ri2r42 

ri2r42 + ri2r42 + r2ir42 

<^33(0)>, 

[l-<T33(0))]ri2r24 

ri2r42 + ri2r42 + r2ir42 

0, / ^ ;, 



(46) 



where (t33(0)) = is the initial population of state \A^}. It is 
obvious that the steady state of the two coupled TLSs depends 
on its initial state. Therefore, the two coupled TLSs exhibit a 
phenomenon of anti-thermalization in the sense that the heat 
bath can not erase totally the initial information of the two 
TLSs. For example, when initially the two coupled TLSs are 
prepared in state \At,}, they will stay in I/I3) forever However, 
in the subspace spanned by the three eigenstates \Ai}, I/I2) , 
and I/I4), the initial information of the two coupled TLSs can 
be totally erased, which can also be seen from Eq. (06]) when 

<T33(0)) = 0. 



E. Steady-state entanglement between the two TLSs 

In the CHB case, in addition to the dipole-dipole interaction 
between the two TLSs, the common bath can also provide a 
physical mechanism to entangle the two TLSs. According to 
Eqs. (I33]) and ( |34] |. the steady-state density matrix elements of 
the two TLSs in the CHB case have the same form as those 
given in Eq. ( |35] |. However, now the steady-state solutions of 
the eigenstate populations (ryy)jj are given by Eq. (|39T l. Cor- 
respondingly, we can obtain the concurrence between the two 
TLSs in terms of Eq. ( |36] |. For a given nonresonant 6, the fig- 
ure in the CHB case is very similar to Fig.|4](so it is not shown 
here). The phenomenon of threshold temperature also exists 
in the CHB case. 



weaker than the TLS-bath couplings, respectively. In the IHB 
case, we have found that, when the two IHBs have the same 
temperatures, the two coupled TLSs could be thermalized in 
eigenstate representation with the same temperature as those 
of the heat baths. However, in the case where the two heat 
baths have different temperatures, just when the energy de- 
tuning between the two TLSs satisfies a special condition, the 
eff'ective four-level system formed by the two coupled TLSs 
can be thermalized with an immediate temperature between 
those of the two heat baths. In bare-state representation, we 
have found a counterintuitive phenomenon that the tempera- 
ture of the TLS connected with the high-temperature heat bath 
is lower than that of the other TLS which is connected with the 
low-temperature heat bath. In the CHB case, the two TLSs in 
eigenstate representation could be thermalized with the same 
temperature as that of the heat bath for nonresonant cases. In 
bare state representation, we have found that the TLS with a 
larger energy separation can be thermalized to a thermal equi- 
librium at a lower temperature. We have also found a phe- 
nomenon of anti-thermalization of the two TLSs in a common 
heat bath in the resonant case. In addition, we have studied 
the steady-state entanglement of the two TLSs in the IHB and 
CHB cases. It has been found that there exist threshold tem- 
peratures for the steady-state entanglement generation. 

Finally, we present some discussions on the thermalization 
time over which the thermalized systems evolve from their 
initial states to steady states. Mathematically, the thermaliza- 
tion time for a system should be infinite because the long-time 
limit (lim,^oo) is needed to make sure that these density ma- 
trix elements evolve to their steady-state values. From the 
viewpoint of physics, we might introduce some time scales 
to describe a thermalization, as the half-life of an exponen- 
tial decay. However, for present systems, the evolutions of 
these density matrix elements [i.e., the transient solution of 
Eq. ( fT0t 1 are not purely exponential functions. At the same 
time, these evolutions depend on the initial conditions. There- 
fore, it is needed to introduce the time scales under given ini- 
tial conditions, other than a universal time scale for a quantum 
thermalization. 



V. CONCLUSION AND DISCUSSIONS 

In conclusion, we have studied the quantum thermaliza- 
tion of two coupled TLSs which are immersed in either two 
IHBs or a CHB. We have characterized the temperatures of 
the two coupled TLSs in eigenstate and bare-state representa- 
tions when the coupling between the two TLSs is stronger and 
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Appendix A: Derivation of quantum master equation ^8) for the IHB case 

In this appendix, we give a detailed derivation of quantum master equation ^, which describes the evolution of the two TLSs 
immersed in two IHBs. In the interaction picture with respect to //o = //jlss + ^^ , the interacting Hamiltonian Q becomes 

H^™\t) = [Ti3Bi3(0 + T24B24(0]e'''" + [TnBnit) + T34B34(0]e'^^' + h.c, (Al) 

where T,y = |/li)(/lj| and we introduce the noise operators 

Bnit) = sin(6i/2)A(f) + cos(6'/2)fi(f), Bi3(0 = cos(0/2)A(f) - sin(6l/2)S(r), 

B24(0 = cos(6i/2)A(0 + sin(6'/2)£(f), ^34(0 = - sin(6'/2)A(0 + cos(6'/2)B(f), (A2) 

with A(t) - YjjSxpf'"'"'' ^nd B{f) - Yjk gikhke^"^'*' ■ Under the Born-Markov approximation (Tl, the quantum master equation 
reads 



dt'TvB [i/f™V), K™V - t'),ps(t) »Pb]] : 



(A3) 



where Tr^ stands for tracing over the degrees of freedom of the baths. We assume that the two baths are in thermal equilibrium 



.(«) «. Jb) 



J") 



7(«)^ 



.(*) 



ih\ 



/(«)^ 



state pB = p;;'' ® p;,7 with pj;;^ = Z;' expi-PiH'^') and pj"' = Z,' exp(-/32//^"0, where Z„ = TrBjexp(-/3ii/^"0] and Z^ = 



^W^ 



TrBjexp(-y62^V)] '^^^ ^^^ partition functions of the two baths, respectively. The parameters ySi - 1/Ti andyS2 = l/r2 are the 

I wave approximation, ■( 



inverse temperatures of the baths for TLSl and TLS2. Through making the rotating wave approximation, we obtain 

: dt'e'^'-^rE^P''{Bl(-t')Bijm-TjiPs dt' cT''^^^-^^'^'' ' 
Jo Jo 



PS ^Y, 

(ij) 



TijPS-i 



+ T 



Xoo r^oo 



{ijM) L 




TijPsTkl 



/-»co /-*oo 

dt'e"^^''^-^'i^''{Bl(-t')Bijm + T,i,psTji A'e'<^''-^^P''<B;.(-f')fi/-t(0)) 
Jo Jo 



+ h.c., (A4) 



where the summation parameter (/, j) in the first line of Eq. ( IA4l i can take (/, j) = ( 1 , 2), ( 1 , 3), (2, 3), and (2, 4), and the summation 
parameter (/;', kl) in the third line of Eq. ( lA4l i can take (ij, kl) - (12, 43), (13, 42), (31, 24), and (43, 12). Here the bath correlation 
functions are defined by {X{t)Y{t')) = TxBV^{t)Y{t')pB], and we use the property {X{t)Y{t')) = {X{t - t')Y{Q)) = {X{Q)Y{t' - i)) 
of the correlation functions. To derive the quantum master equation, we need to calculate the one-side Fourier transform of the 
correlation functions in Eq. (IA4b . For simplicity, in the following we only keep the real parts of the one-side Fourier transforms 
of the correlation functions and neglect their imaginary parts since the imaginary parts only contribute to the Lamb shifts, which 
are neglected in this work. The real parts of the one-side Fourier transform of the correlation functions can be obtained as 
follows: 



and 



Re 
Re 
Re 
Re 

Re 
Re 
Re 
Re 



/-»oo 

£/f'e''="'<B24(0)<(-/')> 
Jo 

dt'e-''^''{Bl^{0)B2A{-t') 
Jo 

Xoo 
dt'e'''''{Bn(0)Bl2(-t')} 

/^OO 

dt'e-'''-''{B\^(Q)Bn{-t')) 
Jo 

dt'e-''^''{B2A{-t')B\^m 
Jo 

dt'e-''^~''{Bn{-t')Bl^m 
dt'e'''''{Bl^{-t')Bnm 

XOO 
dt'e'''-''{B\^{-t')B^m) 



= Re 
= Re 
= Re 
= Re 

= Re 
= Re 
= Re 
= Re 



dt'e"^' {Bu{0)BU-t')} 



= ri. 



y-»CO 

dt'e-'''''{Bl^(0)Bi3(-t')} 
Jo 

dt'e'^^'-'iBumBli-t')} -- 

Jr»co 
dt'e-'''-''{B'l(0)Bu(-t')) 


/-*oo 

dt'e-''^''{Bi,(-t')Bl(Q)) 
Jo 

Xoo 
dt'e-'''-''{B34(-t')B\^m 

dt'e'"''{B\^(-t')B2m) 
Jo 

y-»CO 

dt'e'''-''{Bl(-t')Bnm 
Jo 



= r4 



(A5) 



= Ai, 

= A2, 
= A3, 
= A4, 



(A6) 
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where the parameters F, and A, have been defined in Eq. ([9]). Based on the above one-side Fourier transforms of these correlation 
functions, we can obtain those for other correlation functions. By substituting them into quantum master equation ( IA4b and 
returning to the Schrodinger picture, we can obtain quantum master equation ([8]). 

In the following we give an example for calculation of the one-side Fourier transform of correlation function. 



Re 



I dt' €-•'="'{ 

Jo 



<B ,(0)B24(-f')) 



cos2(0/2)Re 



+ sin^(6'/2)Re 
which is based on the fact that there is no correlation between the two IHBs. We can calculate 



I dt'e-''='''{AHO)A(-t')) 
Jo 

I dt'e-'''"'{BHO)B(-t')) 
Jo 



(A7) 



Re 



dt'e-'"'' 



<A'(0)A(-f')> 



2_^gljna(a>aj)7r6{cJaj - si) = nQ{ea)glisi)naiei) = jaie\)na{si), (A8) 



where we introduce the rate jaisi) = 7Tgaisi)g^{si) and the average thermal excitation na(ei) - l/[exp(ei/ri) - 1]. Note that 



here we have also used the formula 



r°° ' 1 

dt'e^"^' =7r5(w)±/P-, 
Jo w 



(A9) 



where the sign "P" stands for the principal value integral. Similarly, we can obtain Re L dt'e "'''' {B\0)B(-t')} - yh(si)nb(ei) 



with jbisi) - ^Qb(e\)82(si). Therefore, we have 



Re 



dt'e-""' 
Jo 



<B ,(0)B24(-f')> 



= cos^{e/2)ya(si)na{si) + sm^{0/2)yt{si)nb(Bi) = r2. 



(AlO) 



With the same method, the one-side Fourier transform for other correlation functions can also be obtained. 



Appendix B: Derivation of the rates In Eq. f37t for the CHB case 

In the CHB case, the interaction Hamiltonian has the same form as Eq. (lAll i. But now the noise operators become 

Bn(t) = sin(6'/2)Ai(O + cos(0/2)A2(O, Bi3(0 = cos(6i/2)Ai(0 - sin(6»/2)A2(0, 
^24(0 = cos(0/2)Ai(f) + sin(0/2)A2(O, ^34(0 = - sin(0/2)Ai(O + cos(0/2)A2(f), 



(Bl) 



with Ai(f) = Yij8\fij^ "^'' ^nd A2(0 = Yijgiflje '"''. The quantum master equation for the CHB case has the same form as 
Eq. ( IA4l i. but now the real parts of the one-side Fourier transform for the correlation functions become 



and 



Re 

Re 
Re 
Re 



Xco 
e-'''-'' {Bii{-t')B:\^{0))dt' 

e-''"'{B24(-t')Blmdt' 
Jo 

Xco 
e''^-''{Bl2(-t')Bnmdt' 

e'"''{Bl(-t')B24mdt' 





= ri2. 


Re 


= r24. 


Re 


= r2i. 


Re 


= r42. 


Re 



e-''"' {Bu(-t')BUO))dt' 



e-''-^'' {BM{-t')Blmdt' 



f 

Jo 
Jo 

/-*CO 

e''"' {Bl^(-t')Bi3mdt' 
Jo 

e''^-''{Bl(-t')B,4(Q))dt' 





= ri3, 

= r34, 
= r3i, 

= r43, 



Re 
Re 



Re 
Re 



r e-''"''{B24{-t')B\^mdt' 
Jo 

/-*CO 

e-'''-''{Bn{-t')Bl^mdt' 
Jo 

e""' {B\^{-t')Bnmdt' 
Jo 

e''^'-'' {Bl^{-t')Bn{t))dt' 
Jo 



= Re 
= Re 
= Re 
= Re 



e-''"'{Bi3i-t')Blmdt' 
Jo 

e-''^^''<B34(-f')Bl2(0)Mf' 
Jo 

e'^"' {Bl^(-t')B24mdt' 



e'^-'' {B\^{,-t')Bi4{0))dt' 
Jo 



= Ai, 

= A2, 
= A3, 
= A4, 



(B2) 



(B3) 
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where the parameters Tjj and A,- have been defined in Eq. (IJTT i. The one-side Fourier transform of other correlation functions 
can be obtained in terms of the above results. Below, we give an example for calculation of the one-side Fourier transform of 
correlation functions. 



We calculate 



Re 



f 

Jo 



e-^' {BU-t')Bnmdt' 



sin^ (e/2) Re 



+ - sin e^ Re 



e'""'-' {A](-t')Ai(0))dt' 



f 

Jo 

/-*oo 

Jo 



-Hcos^(0/2)Re 



f 

Jo 



e'^'-' {Al(-t')A2mdt' 



{AU-t')A2mdt' 



+ Re 



e''''-' {Al(-t')Aimdt' 



Re 



(B4) 



Xoo "I 

jf'e-2''<A"|(-f')Ai(0)> = ^ 8ljnicoj)7:6(cjj - S2) = 7re(£2)^?(e2)«fe) = 7i(s2)Hs2), (B5) 



where we introduce the rate j\(s2) - ng(s2)gi(s2) and the average thermal excitation n{s2) - l/[exp(s2/T) - 1]. Using the same 
method, we can obtain 



Re 



dt'e"'-'' 
Jo 



{A(-t')A2m 



^ 8ij82jn(ajj)7T6(ajj - S2) = 7Tg(s2)gi(e2)82(e2)n(e2) = ri2(e2)n(e2), (B6) 



where the rate ynCei) = ^gis2)8\(s2)82(s2) = ^J7\is2)72is2)■ Similarly, we have Re [ j^°° dt'e"''-''{A!^{-t')A2iO))j = r2(s2)n(s2) 
and Re [J^ t//'e"'2''<A^(-f')Ai(0))] = ri2(e2)n(e2). Therefore, we obtain 



Re 



f 

Jo 



e''^-''{B'U-t')Bnmdt' 



= fsin(0/2) Vnfe) + cos(0/2) Vr2(e2)]^ nfe) = Tji 



The one-side Fourier transform for other correlation functions ca n also be obtained with the same method. 



(B7) 
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